{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computes and plots geoid undulations from a gravity model\n",
    "\n",
    "#### PYTHON DEPENDENCIES:  \n",
    "- [numpy: Scientific Computing Tools For Python](https://numpy.org)  \n",
    "- [matplotlib: Python 2D plotting library](http://matplotlib.org/)  \n",
    "- [cartopy: Python package designed for geospatial data processing](https://scitools.org.uk/cartopy)  \n",
    "- [lxml: processing XML and HTML in Python](https://pypi.python.org/pypi/lxml)  \n",
    "\n",
    "#### PROGRAM DEPENDENCIES:  \n",
    "- utilities.py: download and management utilities for syncing files\n",
    "- geoid_undulation.py: geoidal undulation at a given latitude and longitude  \n",
    "- read_gravity_model.py: reads the coefficients for a given gravity model file  \n",
    "- calculate_tidal_offset.py: calculates the C20 offset for a tidal system  \n",
    "- real_potential.py: real potential at a latitude and height for gravity model  \n",
    "- norm_potential.py: normal potential of an ellipsoid at a latitude and height  \n",
    "- norm_gravity.py: normal gravity of an ellipsoid at a latitude and height  \n",
    "- ref_ellipsoid.py: Computes parameters for a reference ellipsoid  \n",
    "- gauss_weights.py: Computes Gaussian weights as a function of degree  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.rcParams['axes.linewidth'] = 2.0\n",
    "matplotlib.rcParams['font.family'] = 'sans-serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = ['Helvetica']\n",
    "matplotlib.rcParams['mathtext.default'] = 'regular'\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.ticker as ticker\n",
    "import cartopy.crs as ccrs\n",
    "import ipywidgets as widgets\n",
    "import geoid_toolkit.utilities\n",
    "from geoid_toolkit.read_gravity_model import read_gravity_model\n",
    "from geoid_toolkit.geoid_undulation import geoid_undulation"
   ]
  },
  {
   "source": [
    "#### Choose gravity model"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-- list gfc models from GFZ ICGEM\n",
    "MODELS = geoid_toolkit.utilities.icgem_list()\n",
    "modelDropdown = widgets.Dropdown(\n",
    "    options=sorted(MODELS.keys()),\n",
    "    value='GGM05C',\n",
    "    description='Model:',\n",
    "    disabled=False,\n",
    ")\n",
    "display(modelDropdown)"
   ]
  },
  {
   "source": [
    "#### Read gravity model coefficients"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-- gfc file with spherical harmonic coefficients\n",
    "MODEL = MODELS[modelDropdown.value]\n",
    "GRAVITY = geoid_toolkit.utilities.get_data_path(['data',MODEL[-1]])\n",
    "MD5 = geoid_toolkit.utilities.get_hash(GRAVITY)\n",
    "#-- Download coefficients from GFZ ICGEM server\n",
    "geoid_toolkit.utilities.from_http(['http://icgem.gfz-potsdam.de',*MODEL],\n",
    "    local=GRAVITY, hash=MD5, verbose=True)\n",
    "#-- use maximum degree and order of model\n",
    "LMAX = None\n",
    "#-- use original tide system\n",
    "TIDE = None\n",
    "#-- read gravity model Ylms and change tide if specified\n",
    "Ylms = read_gravity_model(GRAVITY,LMAX,TIDE=TIDE)\n",
    "#-- extract parameters\n",
    "R = np.float(Ylms['radius'])\n",
    "GM = np.float(Ylms['earth_gravity_constant'])\n",
    "LMAX = np.int(Ylms['max_degree']) if not LMAX else LMAX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate map of geoid height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-- PURPOSE: calculate geoid heights at a set of latitudes and longitudes\n",
    "dlon,dlat = (1.0,1.0)\n",
    "lon = np.arange(-180+dlon/2.0,180+dlon/2.0,dlon)\n",
    "lat = np.arange(-90+dlat/2.0,90+dlat/2.0,dlat)\n",
    "nlon = len(lon)\n",
    "nlat = len(lat)\n",
    "#-- reference to WGS84 ellipsoid\n",
    "REFERENCE = 'WGS84'\n",
    "#-- Gaussian Smoothing Radius in km (default is no filtering)\n",
    "#-- no gaussian smoothing\n",
    "GAUSS = 0\n",
    "#-- calculate geoid at coordinates\n",
    "N = np.zeros((nlat,nlon))\n",
    "for i in range(nlat):\n",
    "    N[i,:] = geoid_undulation(np.ones((nlon))*lat[i], lon, REFERENCE,\n",
    "        Ylms['clm'], Ylms['slm'], LMAX, R, GM, GAUSS=GAUSS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create output plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-- setup Plate Carree projection\n",
    "fig, ax1 = plt.subplots(num=1, nrows=1, ncols=1, figsize=(10.375,6.625),\n",
    "    subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "#-- contours\n",
    "PRANGE = (-80,80,20)\n",
    "levels = np.arange(PRANGE[0],PRANGE[1]+PRANGE[2],PRANGE[2])\n",
    "norm = colors.Normalize(vmin=PRANGE[0],vmax=PRANGE[1])\n",
    "\n",
    "#-- plot image with transparency using normalization\n",
    "im = ax1.imshow(N, interpolation='nearest', cmap=cm.viridis_r,\n",
    "    extent=(lon.min(),lon.max(),lat.min(),lat.max()),\n",
    "    norm=norm, alpha=1.0, transform=ccrs.PlateCarree())\n",
    "\n",
    "#-- add generic coastlines\n",
    "ax1.coastlines()\n",
    "\n",
    "#-- draw lat/lon grid lines\n",
    "GRID = [15,15]\n",
    "grid_meridians = np.arange(0,360+GRID[0],GRID[0])\n",
    "grid_parallels = np.arange(-90,90+GRID[1],GRID[1])\n",
    "gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,\n",
    "    linewidth=0.1, color='0.25', linestyle='-')\n",
    "gl.xlocator = ticker.FixedLocator(grid_meridians)\n",
    "gl.ylocator = ticker.FixedLocator(grid_parallels)\n",
    "\n",
    "#-- Add horizontal colorbar and adjust size\n",
    "#-- extend = add extension triangles to upper and lower bounds\n",
    "#-- options: neither, both, min, max\n",
    "#-- pad = distance from main plot axis\n",
    "#-- shrink = percent size of colorbar\n",
    "#-- aspect = lengthXwidth aspect of colorbar\n",
    "cbar = plt.colorbar(im, ax=ax1, extend='both', extendfrac=0.0375,\n",
    "    orientation='horizontal', pad=0.025, shrink=0.90, aspect=22, \n",
    "    drawedges=False)\n",
    "#-- rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "#-- Add label to the colorbar\n",
    "cbar.ax.set_xlabel('Geoidal Undulation', labelpad=10, fontsize=20)\n",
    "cbar.ax.set_ylabel('m', fontsize=20, rotation=0)\n",
    "cbar.ax.yaxis.set_label_coords(1.04, 0.15)\n",
    "#-- Set the tick levels for the colorbar\n",
    "cbar.set_ticks(levels)\n",
    "cbar.set_ticklabels(['{0:d}'.format(ct) for ct in levels])\n",
    "#-- ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=27, labelsize=20,\n",
    "    direction='in')\n",
    "\n",
    "#-- axis = equal\n",
    "ax1.set_aspect('equal', adjustable='box')\n",
    "#-- no ticks on the x and y axes\n",
    "ax1.get_xaxis().set_ticks([])\n",
    "ax1.get_yaxis().set_ticks([])\n",
    "#-- add main title\n",
    "ax1.set_title(modelDropdown.value, fontsize=24)\n",
    "ax1.title.set_y(1.01)\n",
    "\n",
    "#-- stronger linewidth on frame\n",
    "ax1.outline_patch.set_linewidth(2.0)\n",
    "ax1.outline_patch.set_capstyle('projecting')\n",
    "#-- output to file\n",
    "fig.subplots_adjust(left=0.04,right=0.96,bottom=0.05,top=0.96)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}